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Abstract 

We show that computing the coefficients of the Taylor expansion of the solution 
of the off-equilibrium dynamical equations characterizing models with quenched dis- 
order is a very effective way to understand the long time asymptotic behavior. We 
study the p = 3 spherical spin glass model, and we compute the asymptotic energy 
(in the critical region and down to T = 0) and the coefficients of the time decay of 
the energy. 
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It has been strongly stressed in the last period of time that the non-equilibrium dy- 
namics of glassy systems displays very interesting phenomena. One of the most remarkable 
ones is aging [1]. The dependence on time of measurable quantities does not vanish when 
the time goes to oo, and a true equilibrium is never reached. 

Let us give a simple example of an aging behavior. In spin systems the spin-spin 
correlation function can be defined as 

C(t w , t) = N- 1 ]T ai(t w )ai(t w + t) . (1) 

i 

We can consider a system that has been kept at very high T for times smaller than t — 0, 
and has been quenched to the measurement temperature T at time t — 0. We wait a 
time t w after such a quenching, and we measure spin-spin correlations starting at t w . The 
relevant point is that contrary to equilibrating systems, for large t w and t the dependence 
of C on the waiting time t w does not disappear. In the simplest cases it is found that for 
j- ^ 1 C(t w ,t) can be written as 

C(t w ,t) = f(U, (2) 

where / is not a constant function. Even for asymptotically large times the time transla- 
tional invariance of the correlation functions is never recovered. 

Aging is quite a widespread phenomenon in short range models. Aging in disordered 
systems has the interesting peculiarity of being already present in the mean field approx- 
imation. It is obvious that we can have true, complete aging only in an infinite system. 
One should investigate the dynamics of an infinite system. In this case one can write a 
closed set of dynamical equations([2, 3]) which are the appropriate generalization of the 
famous mode-mode coupling equation for glasses. Unfortunately these equations are cou- 
pled integral equations in two variables and the analytic solutions found so far [3, 4, 5, 6], 
although they give a coherent picture of aging, do not allow to study the time dependence 
of the different quantities. It is then highly desirable to have efficient methods to integrate 
the equations numerically. 

An example of these equations for the spherical spin model with a p-spin interaction in 
the case p = 3 is given by a set of non-linear integro-differential equations ([2, 3]), which are 
obtained by transforming eqs. (14,15) back to their constrained form (see the discussion 
in the following). The same equations also appear in the off-equilibrium dynamics of 
the Amit-Roginsky model [7], where although quenched disorder is absent [8] the Mode 
Coupling Approximation gives exact results. 

It is possible to study such equations numerically [4], but it is rather difficult to study 
the solution for large values of the time. The main problem is of the amount of information 
which needs to be stored, since if the discretization step in time is e the number of numbers 
needed to code the correlation function is L 2 where L = t/e. 

The dynamical equations can be derived starting from the usual Langevin equation for 
the p-spin model 

1,7V 

&i(t) = -n(t)<Ti(t) + J i,i2,.-,i P (7 i2( t )- (7 i P ( t ) +Vi(t) , (3) 

t2<...<i v 

where the J couplings are quenched Gaussian random variables, and r\ is a white noise with 
covariance 2T. For all times t the value of //(£) is chosen such to implement the spherical 
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constraint 



f>(t) 2 = iV. (4) 
i=i 

For mainly technical purposes, and some extetical considerations, we find more suitable 
to switch to a formulation where the spins Sj(t) are unconstrained. All what we have 
have done here can be repeated, with essentially the same degree of complexity, in the 
formulation where the spins are the original constrained variables. 
If we consider the transformation 

Si (t) = N(t)a l (t) , (5) 

with 

M(t) = e J>'^') , (6) 

the dynamical equations become 

1,N 

M(ty-%(t)= ]T ^....^(O-^W+^r^W- (7) 

i2<...<i p 

Defining an effective time r from the relation 

dt = U p - 2 dr , (8) 



we get that 



()S ' (T) - E J^,-^(r)-s ip (r)+^(T), (9) 



® T i2<-<i P 

where the ^ are a gaussian noise with zero mean and variance 



( i (r)( J (r')=2TAfW-r%, r (10) 
The unconstrained system of the s variables is kept at a temperature 

T eff (r)^TM^r), (11) 

in order to keep the original, constrained variables, at a fixed, time independent tempera- 
ture T. 

From the previous equations (9,10) after standard manipulations we obtain a set of 
closed dynamical equations for the correlation function 

C(r,rO = ( S ,(r) Si (r')) , (12) 

and for the response function 

^ - S ( s i( T )) /-.ON 

fl(T,T) = W {) 

For t > t' they are respectively 
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§pC(r, t') = l -p{p - 1) jf dr"C^\r, t")R(t, t")C(t", t') 

+ \v f dr"C^\ry)R{r'y) , (14) 
I Jo 

and 

^-R(r, t') = \p(p - 1) jT dr"C p -\r, t")R(t, t")R(t", t') . (15) 
For equal times R — l, and C verifies 

l^C(r,r) = l -p*[ dr"C^(r,r")R(r,r") 

+ TC(t,t) p * , (16) 

where we used the fact that 

JV(T) = y/C(T,T) . (17) 

At the end of the day we have rewritten our constrained system as an unconstrained 
system, kept at a time-dependent temperature T e jj(r). We stress that this transcription 
has been allowed by the fact that our Hamiltonian is homogeneous. We will discuss the 
system in such an unconstrained formulation. We remind again the reader that all our 
results can be easily reobtained in the constrained formalism, where one can also deal with 
non- homogeneous Hamiltonians as the ones in [4]. 

The method we suggest and employ here to study the asymptotic behavior of the 
dynamic equations (14) and (15) is based on computing the coefficients of the Taylor 
expansion of C and R, defined by 

C(t 1 ,t 2 ) = J2°k,j T l T 2 , 

R(ri,r 2 ) = E^M- (18) 

k,j 

We are always assuming T\ > r 2 . Equations (14) and (15) translate in two coupled iterative 
relations for the Taylor coefficients c and r. Knowing the lower order coefficients we can 
determine the higher orders one. The initial conditions are c^o = 1 and r 0j o = 1- We 
determine in turn the coefficients Ckj and r k j with i + j = u for u = 2, 3, 4.... It takes 
a few hours of a RISC workstation to go up to uj = 100 (even if we analyze eventually 
only coefficients going up to the order 48, see later). We repeat the procedure for different 
values of the temperature T. The complexity of the computation increases as uj 3 . 

In this note, as an example, we give the results of the computation of the energy, 
which happen to be in agreement with the theory [2, 3]. In a forthcoming publication we 
will address the problem of verifying more striking issues of the theory as the asymptotic 
violation of time translation invariance and of the fluctuation dissipation theorem. 

We start the analysis from the Taylor expansion for C(r, r). If we define 
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c 



k 

(d) 
k 



E c H,fc. ( 19 ) 



3=0 

we have that 

oo 

C(T,T)=Y,4 d) T l ■ (20) 
1=0 

In the self-explaining formalism we will use in the following we will be working on the 
Taylor coefficients 

{C {d) (r)} k . (21) 

We have computed the coefficients of the series expansions which follow by using a simple 
program written in C-language. We have at first computed the coefficients 

{Af(r)} k = {C«(r)*} fc = { e C )dt '^} k . (22) 
The coefficients of /i as a function of the unconstrained time r are easily computed as 

, . , 1 d/V, r 1 dAfdr^ , d/V -1 , , nn . 

« T »' = = { M If *>» = { — }i ' (23) 

and for the constrained time t as a function of the unconstrained time r 

M^{f^(r)} fc i. (24) 
By inverting and composing these series expansions we eventually obtain 

M*)}k , (25) 

basically the energy as a function of the constrained time t. It is useful to define at last 
the Taylor expansion of the function 

Under the assumption that, for t — > oo, 

/i(t) - /ioo - , (27) 

one has that 



->-(« + 1) . (28) 

Different techniques may be used to extrapolate the function /3(t) to t = oo. After a 
series of tests we have found convenient to use Pade approximants. We have used the 
diagonal Pade approximants of order 48 to check this behavior 1 . We compute a from j3(t), 
with an high enough value of t. having our best estimate for a we compute /i(t) and the 
extrapolated value /j,(t = oo). Once fi(t) is known, the energy E(t) is obtained by the 
simple relation fx(t) = T — 3E(t) [2]. 

1 Lower order Pade approximants have been used to check stability. 
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Figure 1: Theoretical result for the asymptotic dynamic energy E(T) (solid curve) and our 
values (filled dots). The dynamic critical temperature is T D — |. In the high T region 
E{T) = 1 

Let us start with a brief summary of our results. Our method confirms that the asymp- 
totic dynamic energy, E(T) is correctely predicted by the theory. We find in addition the 
remarkable result that the value of the asymptotic dynamic energy of the dynamics at zero 
temperature is the smooth limit of the values for T ^ 0. This does not seem to be true in 
other spin-glass systems as for example the SK model [9]. We also get a very good estimate 
for the exponent a for T < T D = 0.5. These results have been obtained with very small 
computational effort. 

In figure (1) we compare the theoretical result for the asymptotic value of the dynamic 
energy to our findings. Here the critical temperature is obtained by solving the equation 

7^q D (l - qo) 2 = l,q D = q(T D ) = \ (29) 
that gives T D — |, and the dynamic energy is defined as 
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Figure 2: a as a function of T. In these units the dynamic transition temperature T D — \. 
The smooth curve is just used for presenting the data, and does not have any physical 
meaning of a best fit. 



-^(l + g 2 (l-2g)) 



(30) 



The agreement is very good at low T. The convergence (in time) becomes slower when 
approaching the critical point, and exactly at T D we have the higher discrepancy from the 
exact value (of more or less 2%). In the high T phase the convergence to the asymptotic 
result becomes very fast again. 

In figure (2) we plot the exponent a we have estimated from Pade approximants as a 
function of T. We do not have a precise estimate of the errors (this is a known drawback 
of the Pade approach), but judging from the dispersion of different diagonal approximants 
they appear to be of the order of magnitude of a few percent. If we try to estimate a 
for T > To, in the high temperature phase, where we expect the correlations to decay 
exponentially, we find a value that starts to increase with T and dramatically explodes 
for high T. The relatively large discrepancy of the extrapolated energy value close to To 
makes us suspicious about the possible presence of tricky confluent singularities. Obviously 
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Figure 3: Result for the asymptotic dynamic energy E(t) from the direct numerical in- 
tegration of the dynamical equations (solid curve) and our result as a function of time 
(dotted line), here T = ^f-, i.e. T = 0.125 with our normalization. 

such a problem would reflect itself in the value we are estimating for a. 

The fact that our method works very well is also clear from figure (3), where we com- 
pare the results of our expansion to the result from the direct numerical integration of 
the dynamical equations. We have a perfect matching, but, as we discussed, the numer- 
ical integration is limited to small times, and it is difficult to extract from there reliable 
asymptotic values for the coefficients. 

Summing up, we have seen that the numerical integration of the off-equilibrium mode 
coupling equations can effectively be performed computing the coefficients of the Taylor 
expantion of the correlation and the response functions. As an example of application we 
have computed the energy as a function of time and temperature. In the low temperature 
phase (and down to zero temperature) we find power law decay to the value predicted by 
the mean field theory. 
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